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A strategy for finding transition paths connecting two stable basins is presented. 
The starting point is the Hamilton principle of stationary action; we show how 
it can be transformed into a minimum principle through the addition of suitable 
constraints like energy conservation. Methods for improving the quality of the paths 
are presented: for example, the Maupertuis principle can be used for determining 
the transition time of the trajectory and for coming closer to the desired dynamic 
path. A saddle point algorithm (conjugate residual method) is shown to be efficient 
for reaching a "true" solution of the original variational problem. 



2 



I. OVERVIEW 

In many of the complex systems that are encountered in physics, chemistry and biology, 
the potential energy surface exhibits deep minima separated by large barriers. In such cases 
the dynamics is characterized by long periods in which the system stays in one minimum, 
followed by a jump of very short duration into another minimum. These transitions from 
minimum to minimum are rare but crucial since they reflect important changes in the system, 
such as chemical reactions or conformational modifications of molecules in solution. The time 
interval between these rare but important events can easily exceed the limit of molecular 
dynamics (MD) simulations, since it depends exponentially on the barrier height. In the 
case of a smooth potential energy surface (PES) transition state theory [0 can be used 
and a variety of methods have been devised to locate the PES saddle points @. However, 
in a complex system the PES is rough and exhibits a very large number of saddle points, 
thus calling into question the very notion of transition state |l|, Given the relevance 



of this problem, several ways of studying the PES in complex systems have been proposed 
[3], |5| £| [7j|. These methods are based on the idea of accelerating a traditional MD 
through modification of the PES in uninteresting regions, through a suitable rescaling of 
temperatures, or via an increase of atomic masses in hydrogen-rich systems M . 

A different approach has been proposed by Pratt |J. Here the idea is to sample the 
ensemble of paths that connect one minimum with another. This idea has been developed 
into a powerful scheme by Chandler and collaborators The resulting transition path 



sampling (TPS) method has shown great promise in the study of rare events [11, 12, 13 . 

In this approach a dynamical path that connects two minima is determined. Starting 
from this initial path a sound statistical mechanics path sampling procedure is put into 



place that allows transition states to be determined. The procedure for finding the initial 



path is, however, far from simple and a variety of tricks \T1, 0] have been used to this end. 

Motivated by these considerations, we develop here a procedure that is capable of ob- 
taining a realistic dynamical path once the initial and final positions are given. This is at 
variance with other path-based methods which aim at locating the TS or finding the min- 
imal energy path ||14|| . We expect two benefits from our approach: one is to find reaction 
mechanisms in an unbiased way, the other to extract appropriate reaction coordinates in a 
complex system. The basis of our approach is the variational principles of classical mechan- 
ics. In particular, the Hamilton principle states that every classical trajectory that starts 
from a configuration q A and ends in q B after a time r renders the action (which we call 
Sh) stationary. From this principle, in fact, Lagrange equations of motion can be derived 



15| , [Tj] . It has been shown that only for small r or in the trivial case of a free particle is 
the stationary point a minimum [ jT7|| . 

This poses serious problems in numerical applications of the principle. In fact, algorithms 
for locating saddle points are either computationally very demanding or have a very small 



radius of convergence. For instance, Cho, Doll and Freeman fl8fl , and Elber et al. [|19| , |20 
have proposed a method based on the minimization of the squared norm of the gradient 
of the action. However, this approach to saddle points is numerically challenging, since it 
amounts to worsening the condition number of the problem, and it requires the evaluation 
of second order derivatives. 

Another important variational principle that does not require a priori knowledge of r is 
the Maupertuis principle. In this principle the energy is preassigned and an action is defined 
that is stationary relative to the geometrical trajectory. This is complemented by a relation 
between the geometrical trajectory and the time elapsed, which allows r to be determined. 



Our strategy for finding transition paths is the following. We estimate a reasonable value 
of t and start from the Hamilton principle. The saddle point problem is circumvented by 
adding a penalty function which imposes energy conservation. The resulting trajectory is 
then used as an input for the Maupertuis principle and a new and improved value of r is 
extracted. With this new value of r one can go back to the Hamilton principle and repeat 
the procedure until convergence. The resulting trajectories are extremely accurate, so much 
so that one lands in a region where a quadratic expansion of Sh is valid, and numerical 
methods can be applied to the unmodified Hamilton action in order to find the saddle point 



0. 



This paper is organized as follows: in section [I] we introduce our method based on the 
Hamilton and Maupertuis principles; in sections |H] and [Hj we describe in detail the compu- 



tational algorithms we use; in section IV we present a working example: the application to 



alanine dipeptide; finally, section [V] is devoted to our conclusions. 



II. HAMILTON PRINCIPLE AND MAUPERTUIS PRINCIPLE 



In a previous paper [p2|, we introduced an action-based optimization procedure, relying 
on the well-known Hamilton principle. This principle can be stated as follows. Let us 
consider a classical system with 3N degrees of freedom, characterized by its Lagrangean 
L{q, q) = T — V, where T and V are the kinetic and potential energy respectively. If the 
multidimensional Lagrangean coordinates q(t) are fixed to have the values q(0) = q A and 
q( T ) = qB a ^ the initial and final time and r, then the dynamical path can be obtained 
by making stationary the action 

S H = f L(q,q,t)dt (1) 
Jo 



for all the variations that keep the extrema values q A and q B fixed. Indeed, Lagrange's 
equations of motion are derived directly from ([!]). The stationarity does not guarantee that 
Sh is either a minimum or a maximum. On the contrary: only for sufficiently low values 
of r is the stationary point of Sh a minimum ||17|| , and Jacobi ]23| has shown that it can 



never be a maximum; in fact the most common occurrence is a saddle point. We now 
consider a dynamical path q(t) obtained by integration of Lagrange's equations. According 
to Hamilton's principle the Sh is stationary and we can consider the behavior of the Hessian 
S 2 S as a function of t. As discussed above, for small t the Sh is a minimum and therefore 
the Hessian has all positive eigenvalues. As the system evolves this character is lost and the 
time at which the first zero eigenvalue of the Hessian sets in defines the conjugate point. It 
can be proved that, once q(t) has passed through this point, the trajectory is no longer a 
minimum; for larger t, new positive eigenvalues change sign. In the harmonic oscillator, for 
example, a new conjugate point occurs at each half period of oscillation. Only if the system 
has a sufficiently small number of conjugate points will the number of negative eigenvalues 
be much smaller than the dimension of the space, and these directions could be individuated 



and cured, as shown in |[24| . This is not the case for trajectories of interest, as we will show 



in detail in a forthcoming paper [25 



If the system is conservative the total energy E is a constant of motion. In such a case, 
Maupertuis has derived a variational principle in which the time does not appear explicitly 
(historically, the so-called least action principle of Maupertuis came before Hamilton's prin- 
ciple). If we write the geometrical trajectory in the parametric form q(s), where q(0) = q A 
and q(cr) = q B , the Maupertuis action for a system of N particles, with 3N degrees of 



freedom of mass m i} in the so-called Jacobi's form, reads: 



and Sm is made stationary. In (U) the total time r can be calculated as 



Sm is plagued by the same problem as Sh since its Hessian does not have a definite 
character; furthermore (E — V) is not positive definite, which can cause problems in the 
search for a stationary state. In spite of these difficulties, we will show later how to use this 
principle to obtain realistic paths. 

A. Inverted potential and second order methods 

The occurrence of conjugate points during the time evolution of a system represents 
an obstacle to the development of efficient optimization algorithms. Almost 20 years ago 
Gillilan and Wilson |24| (GW) discussed how to overcome this problem in the search for 
dynamical paths. First of all, GW introduced a discretized Hamilton principle (|l|), with the 
integral substituted by a sum, and the velocities defined by: 



so that the discretized action Sh reads 



(5) 
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Applying the condition of stationarity to the function (||) one obtains a set of discrete 
equations of motion, corresponding to the well-known Verlet algorithm: 



q d+D = 2 n) _ q d-D _ ( q m (6) 

GW pointed out that the problem defined by Sh is isomorphic to that of a polymer on 
a surface. The points along the path are the beads of the polymer and the kinetic energy 
term provides the harmonic force that holds the polymer together. Eq. expresses the 
equilibrium between the harmonic forces and the forces coming from the surface. In this 
isomorphism the saddle point nature of the stationary point is reflected in the unstable 
equilibrium between the elastic force and the derivative of the potential V (q^)- Changing 
the sign of V leads instead to a stable situation and Sh has a minimum. By changing the 
sign of V(q) one recovers the elastic band method. 

Several methods for finding transition states and reaction paths |14, |19[ have been de- 



veloped on the basis of this property. Continuing with the polymer analogy, these methods 
often fail to keep the beads evenly spaced along the path; in particular, a consequence of this 
fact is that the polymer does not pass exactly through the transition state, but overestimates 
the saddle point energy by cutting the PES ("corner cutting") around the saddle point. 

From a more formal point of view one can observe the correspondence between the in- 
verted potential path and the instanton semiclassical theory. On this basis one can describe 



corner cutting as a manifestation of quantum tunneling |26|. If one needs to locate the tran- 



sition state, this effect is undesired and the nudged elastic band method has been devised 



to correct for it 14 



A standard method to search for a stationary point is to look for the minimum of the 
sum of the squares of the derivatives. When applied to our case, the object function Sqm 
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to be minimized is: 



*>« - E (5) = E ('. 0+1) - 2 «. m + ^ + („«)) , (7) 

which becomes a minimum when Sh is stationary. Interestingly enough this action coincides 



with the discretized version of the Onsager-Machlup action (OM-action) introduced in [18|, 
I5 | from a totally different point of view. Minimizing this action requires the computation 
of the second derivatives of the potential energy, which can be rather demanding and even 
prohibitive for ab-initio simulations. Furthermore, since the effective condition number of 
the problem is squared, obtaining accurate solutions is difficult and the resulting trajectories 
exhibit poor energy conservation However, since Som measures the deviation of the 
trajectory from a Verlet one, we shall use its value to assess the quality of a path. 



B. Our method 

A consequence of the Hamilton principle is energy conservation along the path (for sys- 
tems where the Hamiltonian does not depend explicitly on time). In a discretized path that 
obeys @ at each time step, energy conservation is no longer exact (and becomes worse and 
worse as A is increased), but will fluctuate about an average value with a certain variance. 
In a first approximation, this effect can be produced by adding to the original action a 
potential that keeps the energy oscillating around an equilibrium value with a given force 
constant. We decided therefore to introduce and minimize the following modified action Sq 
E3: 



S e = l L(q,q,t)dt + fi (T + V-E) 2 dt, (8) 
Jo Jo 
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where only the ratio between \x and 7 matters, and 7 can be positive or negative. In 
the following we restrict the choice to the two values 7 = ±1. The latter term ensures 
the conservation of energy. In the non-discrete case, the second integral is exactly zero for 
physical paths, whereas it will be small and positive if the integral is discretized. Imposing 
energy conservation, although redundant, has the advantage of driving the minimization 
algorithm toward a reasonable zone of the path space in a natural manner. The value of E 
in (|D can be also used as a parameter in the first stages of minimization, in order to guide 
the system toward a value of total energy compatible with the total time r. Within this 
scheme it is also possible to implement constraints involving other conserved quantities like 
total linear momentum and angular momentum; in this paper we will refer solely to energy 
conservation. 



Czerminsky and Elber |27j have already noted in passing that imposing energy conser- 
vation through a particular choice of parameters in their method could be useful but have 
not subsequently used this remark. Here instead energy conservation plays a major role and 
it is the crucial ingredient for arriving at correct dynamical trajectories in an efficient and 
numerically stable way. 

From the mathematical point of view, the consequence of the second term in (§) is that 
the new functional has a minimum for positive values of /x > /i*. We have demonstrated 



this property empirically in ref. [22 for simple cases. Moreover, apart from the condition 
[A > the value of \i is not critical for the location of the solution; this behavior has also 
been observed in more complicated systems. 

A solution obtained from the minimization of (|8]) will in general differ from the solution 
of the Hamilton principle. In other words, one has to refer to quality factors such as the 
Onsager-Machlup action Sqm in order to judge whether the solutions are dynamically sound. 
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In the following, we will define: 

• Stf-path as a solution of the original variational problem ([!]) in the continuum limit; 

• VA-path as a path obtained by integration of Verlet equations of motion, obtained from 
a discretized Hamilton principle. For a small enough A, a V^-path can be considered 
a faithful sampling of a continuous S^-path; 

• 0-path as a solution of minimization of the functional Sq. 

The total time of the path can be estimated using physical or chemical considerations. 
However, this initial guess can be further refined by an alternate use of Sq and Sm- A 
similar, but more expensive approach (since it requires the second derivatives of the energy) 



has been applied recently by Elber et al. to the computation of trajectories |50] with very 
large time steps. Our procedure is instead as follows. 

• First, we fix points q A and q B in configuration space and calculate V(q A ) and V(q B ). 
They do not need to be two exact minima of the potential: two points in the basin of 
attraction of the two states A and B are sufficient in this context. A good procedure 
is to sample these points from traditional MD simulations performed in the basins of 
A and B. 

• We need an initial guess for the path. To start from the linear path connecting q A 
and q B is the simplest choice. From this linear interpolation, we can easily obtain a 
rough approximation of a minimum energy path either by inverting the sign of V in ([5]) 
and minimizing that functional, or by minimizing Sq with r set to a computationally 
convenient value, and a conserved energy Ei ow < V A . The minimum energy path 
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(although not dynamical) gives a crude estimate of the barrier. This allows us to set 
a first choice for E, the total energy along the path, slightly above the barrier. 

• A first value for the total time r can then be obtained from (|3]) applied to the minimum 
energy path treated as a geometric trajectory. Through this equation, we obtain a 
non-uniform time distribution of the intervals along the path. 

• We map the set of geometric intervals on this Sm trajectory onto a set of time intervals 
on the 0-path, and we redistribute the points of the path on an even grid. 

• With the t obtained at the previous iteration, we use Sq full minimization (possibly 
preceded by a slight randomization of the previous path) and find a local minimum. 

• This solution is used as starting point for a few steps of Sm partial minimization. The 
path will be modified and r will be slightly corrected. Partial minimization can be 
successfully replaced by a few cycles of the conjugate residual algorithm for finding 
undefined stationary points described later in this paper. 

• The procedure is repeated until a lower threshold for a quality factor of the path 
(such as Som) is reached. During the iterative scheme the value of E can be adjusted 
automatically, treating it as a slow degree of freedom during the minimization of Sq. 
The value of r will change accordingly during the Sm phase of the optimization. 

The quality of the paths is very high and for many purposes the calculation can be stopped 
here. However, we shall show in the following that the quality can be further improved and 
that the exact variational solution of the discretized Hamilton's principle can be obtained, 
that is, a VA-path can be calculated. 
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C. Fourier components and integration 

The paths can be described numerically in various ways. The simplest one is to use the 
discretization used in eq. @. In most of the applications we have instead used a Fourier 
expansion. Following Cho et aZ.|Tj|, we write the trajectory as: 



9(t) = , A+ (?^9x)( + g aWsin( r^ ) . (9) 

T Z — ' T 

71=1 

In such a way we implicitly take into account the boundary condition q(Q) = q A and 
q(r) = q B . The advantage of this choice is that one can start with few components, thus 
obtaining smooth trajectories which represent a good initial guess and avoiding local spikes 
in the coordinates which lead to unreasonable values of the kinetic energy. Since in this 
representation the positions are defined at any t, the velocities can be exactly calculated as: 



r(0 _ q B -QA , : „,„,, "'A// 



— a w cos . (10) 



n=l 

With this choice, even if the path is discretized, particle velocities are defined at every 
point of the trajectory. Although (|S|) defines the trajectories in a continuum fashion, the 
numerical evaluation of the integral that defines Sq requires the use of a discrete mesh. In 
principle the action integral does not need to be evaluated on a P point mesh, but this 
seems the natural choice, since it is compatible with the path representation. In this way 
we have only one convergence parameter, and the larger the value of P, the better the path 
description. With this choice one has: 

p 

S e = $>>A {jL({q^}, {q«}) + p(T® + V» - Ef) , (11) 
where A = and the weights Wi depend on the integration algorithm. For instance, using 
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Simpson's rule, which will be our default choice, w — 1/2, Wi — 1 for i — 1 . . . (P — 1), and 
w P = 1/2. 

III. MINIMIZATION AND OPTIMIZATION ALGORITHMS. REFINING THE 

TRAJECTORY 

Since Sq can have several minima, in order to minimize our Sq action we use not only a 
conjugate gradient (CG) algorithm, but also a simulated annealing (SA) procedure. In fact 
CG has the disadvantage of leading toward a nearest minimum; with a simulated annealing 
MD, we assign a fictitious mass to the Fourier components and a "temperature" to the 
system, first evolving it at high temperature in order to better explore the parameter space, 
and then "cooling" it toward a minimum of the Sq. This minimum can subsequently be 
refined through the CG algorithm. 

In any case, during the algorithm we have to calculate the derivatives of the action Sq 
with respect to the Fourier components and set them to zero. We write here the explicit 
form for these equations: 




( 



Qb ~ Qa 




p 

Tim 



(m) 



7i I Am Tin Til An 
) — cos 

T T T 



(1 + 2// (T W + V {1) -E)) + 



a 



cos 



T 



F® ({oW}) 



Til An 



(1-2// (T (/) + V (l) - E)) =0 



sm 



T 



(i=1...3N,n = l...P-l), 



(12) 



where F- 1 ; is the i-th component of the force on the ions at time slice /. 
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A. Specific implementations of the algorithm 

We have developed a computer program for performing the minimization of the func- 
tional described in this paper. The program VERGILIUS is an interface connecting 
the evolution of the path to the calculation, at every time slice, and for every step of the 
optimization, of the potential energies V and of the forces F requested by equations (|i~2|). 

At present, VERGILIUS is interfaced to the plane wave MD code CPMD, a quantum 
chemistry code (GAUSSIAN 98 |28||), and to the biomolecular simulation package ORACP9"f. 

The algorithm is naturally parallel, but since we adopt a Fourier representation, there 
are two additional communication phases among the nodes of a parallel machine performing 
a Se-minimization: the distribution of the coordinates from the Fourier components and 
the building of the gradient. To make this phase efficient, we have adopted a fast Fourier 
transform (FFT) algorithm , and the scaling of the overall computational scheme is very 



good. 



B. Quality factors for the solution. Fourier norm 

Once a trajectory is obtained, the question of its accuracy arises. As already stated, if 
the time step is sufficiently small, the closeness of the trajectory to a Verlet trajectory is 
a good criterion: the lower the value of the OM-action, the better the trajectory. In this 
section we will briefly describe this criterion, possible problems arising from it, and some 
procedures for judging the physical soundness of the solution found. 

Since we are optimizing or minimizing our action in Fourier space, another functional 



introduced by Cho, Doll and Freeman [18 can be more appropriate, namely the norm of 
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the gradient of S# (that is, fljjD with /z = 0) in Fourier space (F-norm): 

E(|%) ^(A£«^ raS ^ + f S i,^. (13) 

i,n y^^i / ^ = o i,n I 

Here i is the particle component index (ranging from 1 to 3N) and / is the time slice index 
(ranging from 1 to P — 1). Due to the definition of the velocities, this norm is not exactly zero 
for a VA-trajectory. We have verified, however, that there is a one-to-one correspondence 
between F-norm and OM-action, such that a low value of the F-norm always corresponds 
to a trajectory that is close to a Va path. 



C. Refining: how to go from a B-path to a VA-path. Conjugate residual method 

Once a 6-path is obtained (possibly improved through the Sq-Sm iterative scheme) and 
a reasonably low value for the OM-action found, we can try to reach a stationary point of 
the original variational problem. Having kept the discretization time step sufficiently small, 
that is, comparable with the fastest intrinsic vibrations of the system, we can be confident 
that the VA-path we are aiming at is a good representation of a S#-path. 

For finding the saddle point of Sh, we use an algorithm described in ref. namely 



the Conjugate Residual (CR) method, which is very similar in spirit to the CG, but is 
not limited to solving positive definite linear systems. This algorithm is most efficient for 
sparse systems. The choice of a suitable parameter space is therefore crucial; whereas for a 
5*0 minimization a Fourier component representation has many advantages, we found that 
once the optimization procedure is sufficiently close to a VA-path, returning to Cartesian 
coordinates with velocities defined as finite differences can be wise, since the Hessian of (§) 
is naturally sparse in this representation. We therefore adopt for this refining phase the 
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representation in Cartesian components {q}. 

We assume that we are close to the stationary point and expand the action Sjj around 
the position {q }: 



S H (q) * S H (q°) +J> - g°) ^ 



dq 3 



d 2 S 
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ii.k 



dq h dq k 



(ft-gg)- (14) 



q° 



The condition of zero gradient leads to the linear system 



Ax = b, 



(15) 



where 



.4 



d 2 S 



ii 



ij ■- 



dqidq-j 



OS 



n 



d qj 



x:=q-q . 



q° 



The algorithm scheme is, then, in pseudo-code pl| : 



(16) 



Compute r := b — Axq ,Po '■= r o 

For j = 0, 1, . . . , until convergence Do: 

a 3 '■= ( r j' Ar j)/( A Pj' A Pj) 
X j _|_ • — X j I Ot jP j 

Pj '■= ( r j+l, Ar i+l)/( r j, Ar j) 
Pj+l -■= r J+l + Pjpj 

Compute Apj + i = Arj+i + PjApj 
EndDo 



This algorithm is the counterpart of the conjugate gradient for non-positive definite ma- 
trices. The only matrix-vector product to be performed is the product At. The Hessian of 
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Sh (with the velocities defined as in (|4|)) is a matrix of linear dimension (P — 1) x 3iV (with 
P the number of time slices, and N the number of particles), and has the following block 
diagonal form: 



_ o-b H _ 

A (U,fi),(tj,fj) - (ti) W) - 

I 

A(w tl H tl + w n $1) -Aw tl $1 



... -Awtp^fjl A(w tP H t p + {w tP ^ 1 +w tP )^ I I). 

V J 
Here are time slice indexes, fi,fj are Cartesian components indexes, H ti is the 

Hessian of the potential at the time slice ti and a 3iV x 3N matrix; rh is the vector of the 

masses and Wi are the weights for discrete integration. Despite the sparse nature of this 

matrix, the computation of the matrix blocks is demanding. We chose therefore to calculate 

this product using a first order expansion: 



d 2 S 



H 



dqidqj 
j 



OS 



H 



OS 



H 



q+\q 



(17) 



q-xq y 



3 2A I dq t 

with A small. With this approximation, the algorithm remains sufficiently efficient if the 
solution is not too far from the starting point. Once convergence is reached, a new quadratic 
expansion of Sh can be performed and the procedure is iterated until a desired value of the 
OM-action is found. 

A good preconditioning of A is desirable in order to accelerate convergence. We found 
a good compromise between computational cost and efficiency using so-called "diagonal 
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preconditioning" , which amounts to modifying the original system Ax = b into the system 
M~ x x = M~ x Ab t with the aim of reducing the condition number A max /A m i n of the matrix 
A, and M defined as 

{A hk h = k 
(18) 
h^k 

A convenient way to parallelize the calculation is for instance to assign a time slice to each 
node. Then the diagonal terms of the Hessian can be calculated exactly using eq. (|T7| ) by 
substituting q with a set of 2>N vectors of dimension 3iV x P — 1, with zeros everywhere 
but in P — 1 positions I + 3N(k — 1), k = 1 . . . P — 1. For example, vector q^ has the form: 



/ \ 

1 (3JV- 3) zeros 1 (3iV - 3) zeros ... 



(19) 



V / 

repeated P— 1 times. With this choice every parallel node will return the desired information 
about the diagonal part of A: the application of eq. [17| to the set of vectors leads to 



a set of vectors G^ l \ from which the known off-diagonal term due to kinetic energy must 
be subtracted. At this point the desired matrix Mhk is formed, and the preconditioned 
algorithm can be applied. 

IV. A WORKING EXAMPLE: ALANINE DIPEPTIDE 

The test system we will use in this section is a model of alanine dipeptide, a small 
peptide made of 22 atoms. The PES of this molecule in vacuum shows two main minima: 



an extended configuration, and an axial configuration [ 13| . The barrier in vacuum is about 
7 kcal/mol. The conformational change from an equatorial state to an axial state has been 
studied thoroughly |L3], in particular with a view to estimating free energies both in the 
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vacuum and in solution. As an example, we will adopt here the united atoms (UA) scheme. 
In this model, only two hydrogens are treated explicitly, whereas the others are contracted 
into superparticles with proper mass and charges, according to the OPLS-AMBER scheme 



33| . The resulting UA molecule is made of 12 atoms. We thus have (36-5) independent 
degrees of freedom. 

The scope of this section is to show how a B-path is obtained and how the Sm can be 
used for refining the path. The quality of the trajectory obtained is confirmed by the fact 
that CR algorithm easily reaches a VA-path using the refined G-path as a trial trajectory. 

In order to obtain a realistic path, we used exactly the iterative scheme described in 
Section ||. We started from a linear interpolation between the two minima. We alternated 
S'e-action minimization (with 7 = — 1 and fi = 25000) and partial 5*M-action minimization 
(using both simulated annealing and conjugate gradient algorithms) in order to obtain a 
dynamical path passing through the barrier. 

The results are shown in Fig. [TJ, where we plot the value of the OM-action as a function 
of the iteration steps. The S^f-action refinement allows the value of the quality factor to be 
reduced by one order of magnitude. 

The potential and total energy profiles for this Maupertuis-improved B-path are shown 
in Fig. |[ The total time has been reduced from r = 2ps to r = 1.59 ps. 

A. Refinement 

At this point, we are confident of being sufficiently close to a S#-path to apply the CR 
method and find the stationary path of Hamilton action. 

We take the solution of the combined Sq-Sm strategy (a 0-path with A = 1.9 fs and 
P = 800) and use it as an input for our CR algorithm. 
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The algorithm allows a VA-path to be obtained with a value for the OM-action as low as 
0.5 x 10~ 6 , 10 4 times smaller than the initial value for the G-path. The resulting solution 
can be considered as a true dynamical trajectory, since A is small with respect to the fastest 
oscillation period in this system. The 0-path and the "true" path will diverge increasingly 
in the minimum region of the PES, where the trajectories are chaotic. This reflects a 
fundamental property of reactive trajectories, which exhibit chaotic behavior in the stable 
basins and regularities in the transition region if the energy gap above the saddle is not too 



large p4[ . In terms of our algorithm the chaotic behavior can be explained as follows: near 
a potential energy minimum the forces (and the accelerations) are very small, and a small 
error in the acceleration can cause a change of sign, leading to an incorrect curvature and a 
dramatic deviation of the approximate trajectory from the true one. 

To give a visual measure of the difference between the two trajectories, we plot them on 
a (0 — ij}) PES, where ip and (ft are the two soft dihedral angles of this molecule, as defined 
in ref . [|T^] . It can be seen in Fig. [5] that the two paths pass through the same transition 
state, and are also very similar in the two basins. 

The Euclidean distance between the two trajectories {rl} and {r2}, defined as 



^(^E^EN''-^ (20) 

is in this case equal to 1.15 A/atom/slice. 

Let us now look for a VA-path with A much smaller than that used up to now. Such a path 
can be useful, since in traditional MD the instability threshold for an integration algorithm 
(like Verlet) is well below that for a variational algorithm with two boundaries. Where a 
single transition state is present along the path we can proceed as follows. We consider our 
Va path obtained with the CR, and a small portion around the transition state. Through 
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interpolation, this small portion of the path is refined and A is reduced; CR is then used 
to find the stationary point of this short trajectory. At this point we consider two adjacent 
points in the discrete path. Forward and backward integration using the Verlet algorithm, 
with these points as initial conditions, lead to a reactive trajectory with a very small time 
step. The energy conservation will be of the same quality as in the original O-path. 

Concerning the choice of the parameters in (H), we observed that the choice of 7 = — 1 



leads to more realistic trajectories. We will explain this fact in a forthcoming paper |25 
The value of [i must instead be tuned in order to render the OM-action as small as possible. 
Since there is a wide range of values of /1 (several orders of magnitude) where this quality 
factor does not change appreciably, tuning is not difficult. 



V. CONCLUSIONS AND PERSPECTIVES 

In conclusion, we have built a complete strategy for obtaining reaction paths close to 
solutions of the variational problem of classical mechanics. We have presented an iterative 
procedure based on Maupertuis and Hamilton's principle; as a result, paths are found that 
are sufficiently close to dynamical trajectories to allow a local projection algorithm like the 
conjugate residual method to be applied for extracting the "correct" trajectory from the 
0-path; this procedure is promising also for larger systems. 

The next step in our research will be to pass to phase space, where canonical momenta and 
coordinates have to be considered as independent variables. Discretization of the Hamilton 
principle in phase space and application of Euler equations to the resulting functional leads 
to the velocity Verlet algorithm of molecular dynamics. Moreover, although the number 
of variables is double because of the introduction of canonical momenta, the kinematic 
conditions (first equation of Hamilton) connecting the momenta to the particle velocities can 
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be fruitfully used as further quadratic constraints, leading to another interesting optimization 
problem. 

We wish to thank D. Aktah, A. C. Levi, A. Gusev and A. Laio for fruitful discussions. 
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FIG. 1: Alanine dipeptide: value of the OM-action during the iteration of the Sq-Sm optimization 
scheme. Even iterations correspond to Sq minimization, odd iterations to partial Sm minimization. 



27 




200 400 600 800 1000 1200 1400 1600 

I (ft) 



FIG. 2: Alanine dipeptide: potential energy V and total energy E profile for a path resulting 
from the S&-Sm iterative procedure (refined) described in the text. The total time has been 



FIG. 3: Alanine dipeptide: <f> — tp Ramachandran plot of the Sq-Sm path (dotted line) and of the 
VA-trajectory (continuous line) obtained through the stationary point CR algorithm (see text). 
The two trajectories are very similar, apart from obvious differences especially in the minima 
regions, where the chaotic behavior of the dynamics is more pronounced. 



